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Developing pharmacological strategies for controlling ionizing radiation (IR) -induced cell death is 
important for both mitigating radiation damage and alleviating the side effects of anti-cancer radiotherapy 
manifested in surrounding tissue morbidity. Exposure to IR often triggers the onset of p53-dependent 
apoptotic pathways. Here we build a stochastic model of p53 induced apoptosis comprised of coupled 
modules of nuclear p53 activation, mitochondrial cytochrome c release and cytosolic caspase activation that 
also takes into account cellular heterogeneity. Our simulations show that the strength of p53 transcriptional 
activity and its coupling (or timing with respect) to mitochondrial pore opening are major determinants of 
cell fate: for systems where apoptosis is elicited via a p53 -transcription-independent mechanism, direct 
activation of Bax by p53 becomes critical to IR-induced-damage initiation. We further show that immediate 
administration of PUMA inhibitors following IR exposure effectively suppresses excessive cell death, 
provided that there is a strong caspase/Bid feedback loop; however, the efficacy of the treatment diminishes 
with increasing delay in treatment implementation. In contrast, the combined inhibition of Bid and Bax 
elicits an anti- apoptotic response that is effective over a range of time delays. 

Understanding the mechanism of cellular response to ionizing radiation (IR) damage is important from the 
perspectives of both radiotherapy and mitigation of radiation damage. Cell response to IR involves several 
protein-DNA and protein-protein interactions, as well as the formation of free radicals that alter cellular 
biochemistry \ Cell death usually takes place several hours after radiation injury. Even if the exposure to radiation 
is brief, its effect on cellular biochemistry may be long-lived depending on the strength of IR\ Moreover, several 
proteins that are expressed transiently after radiation damage may trigger downstream responses that are 
manifested long after the original insult. The responses to treatments that aim at alleviating radiation damage 
(or decreasing the susceptibility to apoptosis in damaged cells) depend on the dosage and duration of exposure, 
the treatment timing, and the dynamics of the proteins that regulate apoptotic events. 

The tumor suppressor protein p53 is a main mediator of cell response to genotoxic stress. p53 regulates 
apoptosis via both transcription-dependent and -independent pathways^ ^ in addition to regulating cell/tissue- 
specific response to radiation by apoptosis-independent mechanisms^. The transcription-independent effect of 
p53 is mediated by its translocation to the mitochondria, although the mechanism is still debated (see review^). 
Previous efforts to model cell response to radiation have been in part stimulated by the observed oscillatory 
dynamics, or repeated pulses, of p53 in response to radiation damage^"^. To this end, deterministic methods^"^ 
and, to a significantly lower extent, stochastic simulations^^ have been adopted. Apoptosis itself has been math- 
ematically modeled independent of p53 response to radiation, using deterministic^ ^"^^ as well as probabilistic 
methods^^'^^. Likewise, there have been efforts to establish the link between p53 activities to DNA damage and cell 
fate using deterministic simulations^^'^^ and methods of limited stochasticity^^. 

With accumulating experimental data, we are now in a better position to construct more detailed models for 
p5 3 -mediated signal transduction in response to IR and use them as a platform for evaluating new polypharma- 
cological strategies. Here, we focus on the biochemical network associated with IR-induced apoptosis and 
examine the time- dependence of p53-mediated apoptotic events. Our approach incorporates cell heterogeneity 
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and subcellular localization, and aims at estimating the response to 
targeted therapies following IR. It applies to conditions where certain 
types of molecules are very small in number yet are major determi- 
nants of system behavior. 

We consider several outstanding issues: (i) the significance of the 
oscillatory behavior of p53 in determining the onset of downstream 
apoptotic events in response to IR, (ii) the role of its transcription- 
dependent and -independent activities in regulating cell suscep- 
tibility to apoptosis, (iii) the effect of the interactions involving 
anti- apoptotic Bcl-2 and pro-apoptotic Bax on cell fate, (iv) the 
impact of the positive feedback loop mediated by Bid/caspase-3, 
and (v) the efficacy of various treatment strategies, e.g. how particu- 
lar combination therapies may elicit anti- apoptotic responses to mit- 
igate IR-induced damage. 

Our results indicate that p53 oscillations are insufficient to induce 
apoptosis per se. Activation of Bax on the outer mitochondrial mem- 
brane (OMM) plays a key role in driving apoptosis, and Bid and Bax 
emerge as targets most relevant for regulating apoptosis. In particu- 
lar, truncated Bid (tBid)/caspase-3 feedback loop is a determinant of 
cell fate or treatment efficacy: if the feedback is sufficiently strong, 
inhibition of pro-apoptotic proteins like PUMA mitigates damage; 
but the effect weakens with delay in inhibitor administration. These 
results offer new insights into novel polypharmacological strategies 
for alleviating IR damage. 

Results 

Mathematical modeling of IR-induced apoptosis. We first con- 
structed a mathematical model for the biochemical network 
associated with IR-induced apoptosis. Our model consists of three 
closely interconnected subnets (Figure 1): an upstream p53 module, a 
mitochondrial module and caspase activation events triggered upon 
mitochondrial cytochrome c (cyt c) release. The interactions are 
distributed across three subcellular locations: nucleus (N), cytosol 
(C) and mitochondria (M), indicated by superscripts. 

p53 module. The p53 module consists of the translocation of cyto- 
plasmic p53 (p53^^^) to the nucleus, especially when the cell is 
exposed to stress^^, p53^^^ tetramerization into (p53^^04^°, and the 
transcriptional activation of Mdm2 by p53. The latter involves four 
steps: (i) generation of messenger RNA, mRNAMdmi^^^ (ii) its trans- 
location to the cytoplasm, (iii) translation of mRNAMdm2^^^ into 
Mdm2^^\ and (iv) translocation of Mdm2^^^ to the nucleus, which 
serves as a negative regulator of p53 by promoting its unbinding and 
ubiquitination^\ p53ub^^^ can translocate to the cytoplasm, and 
p53ub^^^ to the mitochondria^^. 

Under severe genotoxic stress, p53^^^ gains transcriptional activity 
(p53a) for pro-apoptotic BH3-only proteins^^ including PUMA and 
mitochondrial outer membrane permeability (MOMP) inducers 
represented by Bax^^. In parallel, it retains its ability to tetramerize 
and activate mRNAMdmi^ and to translocate to the mitochondria. 

The p53-Mdm2 interaction is mediated by upstream DNA- 
damage- signaling proteins such as ATM, Chk2 or Wipl^. DNA - 
double -strand breaks resulting from exposure to IR induce ATM 
autophosphorylation^^, which then disrupt the p53-Mdm2 inter- 
action upon phosphorylating both proteins. We represent ATM- 
associated upstream events as follows: DNA-damage activates 
ATM, which, in turn, stabilizes p53. 

p53-modulated mitochondrial events and release of cytochrome c. 
Mitochondrial events modulated by p53 proceed via two paths: 
transcription of PUMA and Bax, and translocation of p53 to the 
mitochondria followed by binding of p53^^^ to Bcl-2 to inhibit its 
anti -apoptotic action^^. PUMA also binds Bcl-2 on the mitochon- 
drial membrane^^ with a higher affinity than does p53^^^ and can 
displace p53^^^ from its complex with Bcl-2^ thus, freeing up p53 for 
further activity^^. 



Bax translocates between the cytoplasm and the mitochondria^^, 
its retranslocation into the cytosol being mediated by Bc1-Xl^^. Bax^^^ 
interacts with Bcl-2 to form a complex^° that prevents its activation 
(into Bax*) and ensuing oligomerization to form a MOMP pore. 
However, p53^^^ can bind to Bcl-2 stronger than does Bax and dis- 
place Bax, thus countering/alleviating this effect^^ Activation of 
Bax^^^ is facilitated by the locaUzation of tBid to the OMM^\ which 
induces the insertion of Bax into the OMM^^. Bax is also activated by 
p53^^^' and by PUMA'l Bax* oligomerizes on the OMM'^ to form a 
MOMP pore, which, in turn, promotes the release of cyt c^^^ into the 
cytoplasm^ \ Cyt c release is usually considered as the point of no 
return in mitochondria-mediated apoptosis. MOMP pore also 
enables the release of Smac/Diablo^^^ that inactivates the inhibitors 
of apoptosis (XIAPs), further promoting apoptosis^^. 

Events triggered by cyt c release. Cyt c^^^ forms a complex with Apaf-1 
in an ATP -dependent manner, which assembles into the apopto- 
some complex^^ upon heptamerization, and recruits inactive procas- 
pase-9 molecules to activate them into caspase-9 (C9) and catalyze 
the cleavage of procaspase-3 (proC3) to form active C3^^. C3 and C8 
activated by external death signals^^ truncate Bid, resulting in a pos- 
itive feedback loop that amplifies cyt c release^^, while XIAP inhibits 
the apoptosome^^ and promotes the proteasomal degradation of 
C3^°. 

Synthesis, degradation and inhibition of components. Synthesis and 
degradation of monomeric species (not shown in Figure 1; see Tables 
S1-S2) help establish and maintain steady state conditions in the 
absence of stimuli. We consider four drug targets: PUMA, Bid, C3 
and Bax. Their inhibitors are designated as Ipuma» Isid^ Ic3 and lBax» 
respectively. 

Model simulation, calibration and validation. We adopted stoch- 
astic simulations for two reasons. First, the quantity of some proteins 
such as caspase-3 are expected to be extremely low (or non-existent) 
under homeostatic conditions. Second, as Figure 2A shows, 
stochastic simulations reproduce the sustained oscillations of 
p53^^^ and Mdm2^^^ in accord with experiments^'^'^\ while determi- 
nistic simulations result in damped oscillations. 

To establish initial concentrations of proteins that take into 
account of cell-to-cell variability, we first performed a run of 300 h 
to allow the system to reach steady state conditions in the absence of 
IR. Then, we changed the system parameters to account for IR- 
induced perturbations, and allowed the system to evolve. Return to 
unstressed state occurs upon restoring the parameters after At. 
Figure 2B shows the p53^^^ levels as a function of time for transient 
(At = 12 h) and sustained (At = 56 h) exposure to radiation. Note 
that At does not refer to the duration of radiation but to the length of 
time the kinetic parameters are altered from unstressed (normal) to 
radiation (WT) values (see Supplementary Table S3). The transient 
and sustained exposures are analogous to 0.3 Gy and 10 Gy of y- 
irradiation, respectively. In what follows, simulations were per- 
formed under sustained exposure condition. 

We obtained kinetic parameters from experiments whenever 
available and deduced others from the ranges adopted in previous 
computational studies^^'^^. Details on the sources and choices of 
model parameters are presented in the Supplementary Text. These 
were further calibrated to accurately describe previous experimental 
observations^'^'^\ 

Figure 3A shows the comparison of the model-predicted time 
profiles of p53^^^ and Mdm2^^^ (blue lines) with time series western 
blot data^^ (red dots). 1,000 stochastic trajectories were generated 
(under the same initial conditions but with different random seeds) 
and averaged to mimic the heterogeneous cell population behavior. 
Good agreement is achieved between predicted time profiles and 
experimental data. 



SCIENTIFIC REPORTS | 4 : 6245 | DOI: 1 0.1 038/srep06245 



2 



IR 



p53(c) 



DNA damage 



I — p53(N) — 

I - p53a ^ (d53.). Li— ► 



I 

Mdm2(c) M— 

Mdm2(N)--, 
17 / 



Cytoplasm 



I 
I 

t 20 



(P53a74 

13 



■(P53'^')4 Nucleus 

3 



Bax 

Bax(M) 



PUMA 

I 1^23 



Bcl-2 



BcI-Xl(m) 



BcI-XlJ 



^.Bcl-^<^^3.Bcl-2 



PUMA.Bcl-^ 
p53(M) 



Mitochondrion 
33 _^cytdM) 
Smad^)^ 



Bax 



Bax. 



C8 



I 28 

tBid(M) 

I" 

tBId 



Bid 



Smac.XlAP ^^--^ 
apop.C9.XIAP 



"o 



Smadob O Apaf 
cytc.Apaf 

-1; 



XIAP - 



apop.C9 



C3ub 

' 

C3 



^36 

apop 

proC9 

proC3 



C9ub 



Apoptosis 



Figure 1 | Simplified reaction network diagram of the mathematical model. The diagram highlights the major reactions in the model. Basal protein 
synthesis and degradation reactions are included in the model but not shown. The full list of components, reactions and kinetic equations, and parameters 
are presented in the Supplementary Tables SI and S2. Complexes are denoted by the names of their components, separated by a dot. Single-headed solid 
arrows characterize irreversible reactions and double-headed arrows, reversible reactions. Dotted arrows represent enzymatic reactions. The reactions 
computed by sensitivity analysis (Supplementary Figure S2) to play a significant role are shown by red arrows. Among them, the kinetic steps 20 and 10 (or 
associated rate constants and kio) lead to the respective transcription-dependent and -independent activities of p53. 



We next validated the model using single-cell-based experimental 
data^'^^. We compute the histogram of the period and amplitude 
values of Mdm2^^^ oscillations using the 1,000 trajectories simulated 
under IR. As shown in Figure 3B and C, the simulation results match 
the experimentally detected period^ and amplitude^^ of oscillations, 
indicating that the model captures well cell-to-cell variability property. 

Our model also reproduced key observations downstream the 
signaling cascade (Figure 3D). Note that two stochastic simulations 
may lead to different dynamics even with identical initial conditions, 
and conclusions drawn from a single simulation may not be reliable. 
We applied a statistical model checking (SMC) technique^^ in order 
to make predictions with given confidence level. Briefly, we encode 
dynamical properties of the system using a bounded linear temporal 
logic and devised a SMC procedure to check if the system that 



satisfies a given property passes a statistical test with predefined 
type-I (false positive) and type-II (false negative) errors. A detailed 
description can be found in the Methods. Figure 3D summarizes the 
dynamical properties predicted with high confidence. These prop- 
erties are consistent with existing experimental observa- 
tions^'^^'^^'^^"^^. This analysis further validates the model and 
supports its use for predictive purposes, presented next. 

Sensitivity analysis. In order to evaluate the sensitivity of the model 
to initial concentrations and kinetic parameters we performed a 
sensitivity analysis based on the level of caspase-3 activation. The 
results show that the system is robust to perturbations in the initial 
concentrations of major species, in support of the approach adopted 
for defining initial distributions. 
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Figure 2 | Simulation of p53 and Mdm2 dynamics. (A) Comparison of stochastic (left) and deterministic (right) simulations. The time profiles of p53^^^ 
and Mdm2^^^ were simulated using the stochastic approach (left panel) and a deterministic approach with the same kinetic parameters (right panel). 
The stochastic simulation shows sustained oscillations while the deterministic simulation results in damped oscillations. (B) Radiation exposure is 
initiated at t = 0 (upon alteration of kinetic parameters, which applies for a duration of 12 h (left) or 56 h (right))y and leads to p53 oscillations. The 
amount of the radiation doses for 12 h and 56 h cases are comparable to 10 seconds and 5.6 min y-irradiation (^°Co, 1.8 Gy min"^) treatments, 
respectively^. 



The calculated sensitivities to model parameters are presented in 
Supplementary Figure S2. Strong control over the system response is 
distributed among the p53-Mdm2 interactions (k45), activation of 
the nuclear p53 (kn) and ensuing transcriptional activation of 
PUMA (A;2o), translocation/retranslocation of Bax between the mito- 
chondria and the cytosol (k24~), Bax activation (k3o) and caspase-3 
regulation (k^g^ and k42). The reactions with high influences are 
colored red in Figure 1. 

The analysis highlights the significant role of p53-mediated tran- 
scription-dependent pathway: p53^^^ p53a PUMA/Bax ^ 
g^(M) _^ Bax*. The p5 3 -transcription-independent events (invol- 
ving Mdm2 and p53^^0 have moderate effects and become influ- 
ential when they feed into the Bax activation pathway. These 
results suggest that the strength of coupling of the nuclear p53 mod- 
ule to the mitochondrial apoptotic machinery will play an important 
role in determining the cell fate. In the next sections, we focus on the 
effects of the reactions identified here to be most influential. We 
begin with the examination of the relative contributions of the 
p53 -transcription- dependent and p53-translocation-dependent 
activities of p53 to eliciting mitochondrial apoptotic responses. 

Coupling of the p53 module to mitochondrial pore opening. p53a 
and p53^^^ initiate the p53-transcription-dependent and -indepen- 
dent apoptotic events, respectively, by either upregulating pro- 
apoptotic proteins (represented by PUMA and Bax), or interacting 
with antiapoptotic Bcl-2 family proteins. These may lead to MOMP 
opening, cyt c release and caspase activation depending on the 



coupling between the two pathways, which in turn depend on the 
relative rates of (i) transcriptional activation of pro- apoptotic 
proteins by p53^^^ (kn) and (ii) translocation of p53^*^^ to the 
mitochondria (kio). We varied kn and kio, and simulated the time 
profiles of caspase-3 activation (commonly used as an indicator of 
apoptosis) toward elucidating the relative sensitivity of caspase-3 
activation to these two p53-mediated mechanisms - higher kn and 
kio values representing the dominance of transcription- dependent 
and -independent activities of p53, respectively. The results are 
presented in Figure 4 A and 4B. 

Figure 4 shows that the onset of apoptosis is elicited either by 
enhancing the transcriptional activity of p53 (increasing kn) (A), 
or accelerating the p53^^^ translocation to the mitochondria (increas- 
ing kio) (B), while low efficiency of either process is not sufficient to 
cause apoptosis. To make an assessment of the efficacy of either 
pathway on determining the cell fate, we further predicted the res- 
ponse curves of caspase-3 activation to kn and k^. As shown in 
Figure 4C, the apoptotic response to kn is clearly greater than that 
to kio indicating that the p5 3 -transcription- dependent pathway is 
playing a dominant role in determining cell fate. 

Next, we examined the role of p53 oscillation in mediating apop- 
totic response. The p53 level oscillate with a —5.5 h period due to 
radiation exposure (Figure 2B). It is unclear if the apoptotic response 
is mainly determined, or predominantly affected, by the p53 oscil- 
lation frequency. To address this question, we perturbed the para- 
meters (Supplementary Table S3) to produce a condition, under 
which the p53 level oscillates with a lower frequency (LF) (—11 h 
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Figure 3 | Model predictions and validation. The time profiles of p53^^^ and Mdm2^^^ (A), the histograms of the period (B) and amplitude (C) of 
oscillations for Mdm2^^^ under IR are simulated and compared against previous experimental observations. Blue solid lines and bars depict the simulation 
results and red dots and bars indicate experimental data. The data in bluelines and red dots were normalized so that their maximum value was equal to 1. 
The experimental data in (A), (B) and (C) were extracted from Lahav et al, 2004^^ Geva-Zatorsky et al, 2006^ and Geva-Zatorsky et al, 2010^^, respectively. 
(D) The dynamical properties predicted by the model are supported by experimental evidence^'^^'^^'^'^"^^ (type I error = 0.05, type II error = 0.05, see 
further details in Supplementary Material). 
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Figure 4 | Simulated apoptotic response with different strengths of coupUng between p53 dynamics and mitochondrial events. (A) Caspase-3 time 
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denoted as WT. The associated parameters can be found in Supplementary Table S3. 



period) in response to IR (Supplementary Figure SI). Under such 
conditions, the apoptotic response previously observed forku = 2 X 
10"^ s"^ remained unaffected (Figure 4D), suggesting that slowing 
down p53 oscillations does not dampen the p53-transcription- 
dependent apoptotic response. Interestingly, an apoptotic response 
could be triggered even with a moderate transcriptional activity of 
p53 {ku = 10"^ s~\ Figure 4D), which would fall short of inducing 
apoptosis under WT conditions (see panel A). This might be 
explained by the increase of p53^*^^ level (Supplementary Figure S3) 
due to LF condition, which, in turn enhances the p53-transcription- 
independent pathway. These results imply a dependency of the p5 3- 
mediated apoptosis on the p53 oscillation frequency, and also suggest 
that p53 -transcription-independent pathway might be capable of 
defining cell fate when p53-transcription-dependent pathway is 
suppressed. 

A crucial insight that emerges from these findings is that IR- 
induced cell death depends on the coupling between the regulatory 
p53 module and downstream apoptotic cascade of events, and 
among the two regulatory activities of the p53 module, transcription 
of pro-apoptotic proteins appears to be more effective that those 
driven by the translocation of p53 to the mitochondria. 

Next, we examine the response of proteins downstream of the p53 
module toward shedding further light on the significance of particu- 
lar interactions/reactions in determining cell fate. 



Bax activation. The above analysis highlights the importance of 
transcriptional regulation by p53, while also drawing attention to 
the complementary role of transcription- independent mechanisms. 
Bax is actually affected by both mechanisms. The sensitivity analysis 
above indeed highlighted the importance of Bax activation (kso, ^sh 
k74) for initiating apoptosis (Supplementary Figure S2). As shown in 
Figure 1, in the p53-transcription-dependent pathway, p53a 
prompts Bax activation by upregulating the expressions of Bax and 
PUMA which activates Bax either directly or indirectly (via binding 
Bcl-2 and thereby preventing its anti- apoptotic effect on Bax). 
Similarly, p53^^^ in the p53-transcription-independent pathway 
also activates Bax either directly or indirectly via binding Bcl-2. 
Thus, the affinity of Bax^^^ for Bcl-2, accounted for by the ratio 
^25~l^25^^ controls the efficacy of the indirect Bax activation (via 
either p53-transcription-dependent or independent pathways). 

^25~l^25^ widely varies with experimental conditions and cell 
^^^26,27 Values spanning several orders of magnitude, from micro- 
molar to nanomolar, have been reported for this dissociation con- 
stant by different groups (see the Supplementary Text). Figure 5 
shows caspase-3 and Bax^^\Bcl-2 complex concentrations for two 
different k25~/k25^ values. Comparison of Figure 5A and C shows 
that the caspase-3 levels (blue curve, right ordinate) are unaffected by 
variation in the dissociation constants over two orders of magnitude. 
In contrast, vastly differing amounts of Bax^^\ Bcl-2 complex maybe 
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Figure 5 | Robustness of caspase-3 activity with respect to association strength of Bax^^^ and Bcl-2 in the mitochondria. Panels A and C show onset of 
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formed {black curve, left ordinate) consistent with the change in the 
dissociation constant. Similarly, the different dissociation constants 
had an insignificant effect under anti-apoptotic conditions (panels 
5B and D). 

The model is therefore robust to variations in Bax^^^Bcl-2 asso- 
ciation strength. This suggests that cells exhibiting significantly 
different propensities for Bax^^\Bcl-2 complex formation or dissoci- 
ation can still exhibit the same cellular fate. It also implies that the 
indirect activation of Bax (via its dissociation from the complex 
Bax^^\Bcl-2) has only a secondary effect on the overall cell response. 
This conclusion is further supported by results presented in 
Supplementary Figure S4, which shows that knocking down the 
interactions between PUMA and Bcl-2 does not significantly 
decrease caspase-3 levels (Supplementary Figure S4). 

Caspase/tBid positive feedback loop. Bax activation leads to the 
release of cyt c and causes caspase-3 activation. Caspase-3 cleaves 
Bid and the resulting tBid^^^ further activates Bax. In this section, we 
examine the role of this positive feedback loop, Bax* ^ C3 ^ tBid 
Bax*. 

Figure 6 displays the time evolution of PUMA, tBid^^^ and 
caspase-3, after exposure to IR. PUMA {black) mirrors the p53 
oscillations, presumably due to its direct transcriptional regulation 
by p53. Interestingly, PUMA levels decay after termination of IR 
exposure, but there is sustained tBid fluctuations and newly eli- 
cited caspase-3 production even after the decay of PUMA. 
Sustained activation becomes more prominent with increasing 
coupling to p53 transcriptional machinery (high kn. Figure 6B). 
Caspase-3 activation depends on whether sufficient tBid is acti- 



vated by the time PUMA decays. The increase in [tBid] (and 
associated positive feedback to Bax) is critically important for 
sustained caspase activity. Otherwise, PUMA and Bax upregula- 
tion may fall short of triggering efficient cyt c release. Under these 
circumstances, the extent of Bax activation by p53^^^ {kiQ, k29) or 
PUMA {k^o) ni^y become critical. 

The observed importance of the caspase/tBid feedback loop in 
modulating apoptosis warrants further study. Two of the possible 
ways that this strength can differ in cells is through truncation of Bid 
by caspase-3 and through activation of Bax by tBid (controlled by 
^is)- We varied k28 to study the effect of this feedback on cell fate 
(high kii. Figure 6C-D). Unless k28 is sufficiently high {k28 = 
0.1 |iM"^s"^), no activation of Bax^^^ by tBid takes place, i.e., the 
positive feedback loop that promotes apoptosis is decoupled from 
the apoptotic machinery, and Bid has no influence on apoptosis. In 
the intermediate values of k28 = 0.3 |iM"^s"\ caspase-3 reaches a 
level comparable to that attained with k28 = 0.5 |iM"^s"\ although 
the response is slower. Therefore, in contrast to the strength of Bax/ 
Bcl-2 association, the intensity of the positive feedback loop signifi- 
cantly affects the time evolution of [C3] , suggesting that targeting the 
caspase- 3/tBid loop might be an effective therapeutic strategy for 
mitigating radiation damage. 

Efficacies of polypharmacological strategies. In the sections above, 
we have analyzed the dynamics of apoptotic mediators in response to 
IR. We now focus on identifying potential targets and polypharma- 
cological strategies suitable for mitigating radiation damage. Of 
special interest are treatment strategies that are effective even if not 
administered immediately after IR exposure. 
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Figure 6 | Role of tBid activation for sustained caspase-3 activity and effect of the strength of caspase/Bid positive feedback loop in mediating apoptosis 
and its inhibition. (A-B) Time evolution of the concentrations of PUMA, caspase-3, and tBid at in response to long time IR exposure predicted for two 
different coupling strengths of p53 transcriptional activation to mitochondrial pore opening events (represented by kn). PUMA exhibits oscillations 
echoing the behavior of p53. tBid production remains limited {red curve) on the left panel A (low kn) such that caspase-3 levels are negligibly small. 
Enhanced coupling {right panel) yields a sustained increase in tBid levels, which in turn induces a significant increase in caspase-3 levels, via a positive 
feedback loop. Note the ordinate scale difference between the two panels. (C) The caspase-3 time profiles were simulated for three /c28> which correspond 
to different strengths of caspase/Bid feedback loop. Upon decreasing k28y caspase-3 concentration progressively decreases. (D) Feedback loop strength- 
caspase activation response curves. 



Figure 7 illustrates the efficacy of simulated drug treatments admi- 
nistered at two different times after exposure to IR: immediately 
(with a time lag of 15 minutes succeeding IR exposure, Figure 7 A) 
and after a delay of 12 hours (Figure 7B). The treatments involve 
administration of individual inhibitors for four targets, PUMA, Bid, 
caspase-3 and Bax, or combined therapies targeting pairs of these 
proteins. All (virtual) inhibitors are assumed to have nanomolar 
binding affinity to their target proteins. In the absence of treatment, 
the given radiation dose leads to sustained caspase-3 activity. 

Individual inhibitions of PUMA, Bid or Bax are very effective in 
mitigating radiation damage, provided that treatment is available 
immediately after radiation damage. On the other hand, if there is 
a longer delay before the treatment is administered (Figure 7B), indi- 
vidual inhibitions of Bax and Bid still help in abrogating radiation 
damage but PUMA inhibition is no longer an effective treatment. 
The ineffectiveness of PUMA inhibition is due to the fact that a 12- 
hour PUMA activity is sufficient for the activation of the tBid/ 
Caspase-3 positive feedback loop and activation of Bax via tBid 
dominates the apoptotic response at this late stage. 

In contrast to inhibiting PUMA, Bid and Bax, inhibition of caspase- 
3 itself does not prevent apoptosis irrespective of the timing of inhibi- 
tion. This is an intriguing observation since complete abrogation of 



caspase-3 would be expected to inhibit apoptosis (although not cas- 
pase- independent cell death). However, we have a continual produc- 
tion of procaspase-3, and the activation of caspase-3 upon cleavage of 
procaspase-3, which is enabled by Bid and PUMA, dominates the 
outcome. This observation suggests that targeting caspase-3, itself, 
may not be an efficient strategy for preventing apoptotic response. 

All the combination therapies considered here abrogate apoptosis 
completely (e.g., inhibition of PUMA + Bax, Bid + Bax, PUMA + 
Bid), or significantly (Bid + C3) when administered immediately 
after radiation exposure. This might be expected as Ipumaj Isid ^nd 
Ifiax are already effective alone and their joint administration appar- 
ently does not have an unexpected effect. When the inhibitors are 
administered after 12 hours, the combination Ipuma + Ic3 is inef- 
fective, consistent with the failure of these two inhibitors to prevent 
caspase-3 accumulation individually (Figure 7B). However, the time 
dependence of [C3] is interesting: it shows a short-term response 
suggestive of mitigation of apoptosis, followed by an apoptotic res- 
ponse after a day or two. This combination therapy is therefore not 
effective - again, due to the attainment of sufficiently high level of 
tBid to sustain the pro-apoptotic positive feedback loop that is 
enhancing Bax activation. We then simulated the attenuation of 
the strength of this feedback loop by reducing (from 0.5 to 
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Figure 7 | Potential mitigation of radiation damage via individual and combination therapies. Panels (A) and (B) show the effects of therapies based on 
single targets, PUMA, Bid, caspase-3, Bax, or their combination when inhibitors are administered at 15 min (A) and 12 h (B) after radiation exposure. Bid 
and Bax appear to be the most effective targets even when inhibited with a time delay. The control (no inhibition) is shown in both panels. The efficacy of 
PUMA inhibition strongly depends on the timing of the drug treatment. (C) Time evolution of caspase-3 in the absence of inhibitors (black) and in the 
presence inhibitors of PUMA and C3 {pink), and PUMA and Bax {blue). (D) A reduction in caspase/tBid feedback strength {k28) (from 0.5 to 
0.3 |aM~^s~^) rescues the efficacy of PUMA inhibitor administrated after a delay of 12 h. 



0.3 |iM"^s"^). The results confirmed that Ipuma administered after 
12 hours would be effective in mitigating radiation-induced apopto- 
sis provided that this positive feedback loop is weakened {blue curve, 
Figure 7D). 

These results highlight the importance of examining the long-term 
behavior before making an assessment on the efficacy of a treatment, 
the role of timing in ensuring an effective treatment, and the signifi- 
cance of identifying most susceptible targets for efficacy of the 
treatment. 

Discussion 

Recently, Batchelor et at reported that stimulus -dependent temporal 
dynamics of p53 is an important determinant of cell regulation^^. A 
further study suggests different p53 dynamics (repeated pulses vs 
sustained response) may alter cell fate^^. The present study focused 
on a detailed examination of the dynamics of p5 3 -mediated and 
mitochondrial- dependent modules for regulating apoptosis in res- 
ponse to IR damage. It also highlights the significance of the timing 
between genotoxic stress and therapeutic intervention in the context 
of the cell stochastics. 

Our model analysis showed that the p53 oscillatory behavior in 
response to IR is not, per se, sufficient to explain a cell's susceptibility 
to apoptosis; rather, the dynamic coupling between the p53 module 
and mitochondrial machinery is important. Among transcription- 
dependent and -independent roles of p53, the former plays a major 



role in driving apoptosis via transcriptional activation of Bax. The 
present study also sheds light to the strength of caspase/tBid feedback 
loop as a determinant of apoptotic response as well as treatment 
efficacy. 

Several experimentally testable results emerge from the analysis: 
(i) oscillatory behavior of PUMA in response to radiation damage 
(Figure 6A), (ii) elevated tBid activity even after the decay of PUMA 
when increased caspase-3 levels occur (Figure 6B), (iii) inhibition of 
Bid and/or Bax as more effective strategies for mitigating radiation 
damage, compared to inhibition of PUMA (Figure 7), and (iv) the 
need to have more potent inhibitors of caspase-3, compared to those 
of PUMA, Bid or Bax, in order to effectively mitigate radiation- 
induced cell death. 

The last two points are especially relevant to developing efficacious 
therapies to regulate apoptosis. Our earlier experiments and phar- 
macophore modeling studies showed that inhibition of PUMA can 
be an effective strategy for mitigating radiation-induced cell death^^. 
However, in the presence of a strong caspase/tBid positive feedback 
mechanism that amplifies Bax activation, even a potent inhibitor of 
PUMA would fail if not administered immediately after IR-exposure. 
In contrast, for cells with a reduced strength of caspase/tBid feed- 
back, PUMA inhibition could be effective even if administered after 
substantial delay (Figure 7D). The above considerations do not 
account for the structural promiscuity of these targets: PUMA is a 
BH3 domain-only protein and Bax also contains a BH3 domain. An 
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inhibitor of PUMA may also bind and inhibit Bax, thus resulting in 
the effective therapy of combined Ipuma + ^Bax (Figure 7B). 

It is possible to test the hypothesis (iii) above using the following 
two approaches. Firstly, a knockdown of PUMA when performed 
several hours after radiation damage should be ineffective in mitigat- 
ing apoptosis (although when administered immediately, it would) 
for cells with a strong caspase/Bid feedback. Secondly, knocking 
down Bid mRNA several hours after radiation damage should be 
more effective in mitigating apoptosis at both shorter and longer 
time delays. Inhibitors for Bax^°, Bid^^ and PUMA^^ have been syn- 
thesized, that could be used for these testing purposes. 

We have also performed an extensive sensitivity analysis of the 
network presented in Figure 1, and identified robust reactions whose 
parameter changes do not significantly affect the cell behavior. For 
instance, the binding/dissociation constant between pro-apoptotic 
Bax and anti-apoptotic Bcl-2 (A;25"), on the mitochondria, has a 
relatively small effect (Supplementary Figure S2). On the other hand, 
our study highlights the significance of the translocation of Bax from 
the cytosol to the mitochondria (succeeding the transcriptional 
activation of Bax and PUMA by p53) for eliciting apoptotic response. 
The translocation Bax^^^ ^ Bax^^^ {^24') unambiguously emerged 
from our sensitivity analysis to be a key initiator of the highly influ- 
ential interactions that lead to apoptosis, indicated by red arrows in 
Figure 1. This finding is in support of the recently proposed retro- 
translocation of Bax from the mitochondria to the cytosol, mediated 
by Bc1-Xl, as an effective mechanism for regulating pro-apoptotic 
activity^^'^^. Enhancing retrotranslocation may indeed be an essential 
strategy for decreasing susceptibility to cell death induced by IR. 

An important concept that has recently emerged in cancer therapy 
is mitochondrial priming - whereby cells that are exposed to higher 
concentrations of BH3-only proteins such as PUMA, NOXA and 
Bim are designated as Veil-primed'^^. Well-primed cancer cells are 
closer to the apoptotic threshold than are normal cells, leading to 
more efficient responses to therapies that induce apoptosis. 
Although, mitochondrial priming is used in conjunction with 
chemotherapy, the results in Figure 4 show some parallels to the 
priming process. Increased kn beyond a threshold value (panel D) 
presumably leads to overexpression of BH3 -domain-only proteins 
(e.g. PUMA), promoting apoptosis. The same panel further shows 
that a switch from anti- to pro-apoptotic state (indicated by an 
abrupt change in [C3]) is alternatively elicited when the rate for 
the translocation p53^^^ p53^^^ exceeds a threshold value {kio > 
2 X 10"^s"^), indicating that transcription-independent activity of 
p53 alone may also induce apoptosis. Inhibition of p53 binding to 
mitochondria has indeed been reported to protect mice from gamma 
radiation, consistent with the predicted behavior^^. 

The development and application of systems models (whether 
stochastic or deterministic) to understanding cellular fate is an evol- 
ving field. We note that the sustained oscillation of p53 and Mdm2 
can be maintained under deterministic conditions using certain 
parameter settings^^. However, it is often difficult for deterministic 
models to reproduce the observed statistical features (e.g. distri- 
bution of oscillation periods and amplitudes) of a cell population 
(Figure 3B and C), because these quantities tend to be constant under 
deterministic condition. The current model and parameters under 
stochastic conditions capture the cell-to-cell variability observed in 
single-cell experiments^^'*^ while the damped oscillation behavior 
under deterministic condition is also consistent with the work of 
Alon and coworkers^^. Further experimental data on changes in par- 
ticular protein levels in response to radiation exposure will help 
refine the models, as well as establish rate parameters (via direct 
implementation of experimental data, or statistical inferences from 
expression/activity data). In the present setting, we included in our 
model the key events relevant to p53-mediated apoptotic response to 
IR. The basal activity and effects of other mechanisms such as mito- 
chondrial fusion/fission events^^'^^ are implicitly captured by the 



kinetic parameters. A more comprehensive study of the role of p53 
dynamics would require to expand the model to include p53-regu- 
lated cell cycle arrest and DNA repair events as well mitochondria- 
targeted inhibitors of cyt c peroxidase^^ and caspase-independent 
apoptotic interactions^^'^^. These further extensions may help design 
more efficacious polypharmacological strategies for controlling cell 
susceptibility to apoptosis under different disease states and envir- 
onmental challenges. 

Methods 

Extended Gillespie algorithm. We adopt an extension of Direct Reaction version of 
the Gillespie algorithm^", which takes account of the stochasticity of interactions and 
heterogeneity of the cellular environment^^ In this approach, one reaction is chosen 
at a time, and the time is advanced based on the overall reaction propensity at that 
time. Consider the following schematic reaction, 

X + Y ^XY 



The stochastic rate, c, of the reaction is related to the macroscopic kinetic rate 
constant, /c, as c = AiV, where V is the reaction volume; and the instantaneous 
propensity for this reaction is*^°: 

aa = cNxNY = kNxNYV 

where Nx and ATy are the number of molecules X and Y in the reaction volume and 
the subscript a denotes that the reaction index in the system of M chemical reactions 
(1 < a < M). The underlying assumption is that the system is in thermal equilibrium 
(the number of reactive collisions are much less than that of non-reactive collisions) 
and well-mixed. However, cellular systems are spatially heterogeneous: for example, 
the mitochondria and the cytoplasm provide vastly different environments; and, 
there may be heterogeneities within these localizations themselves. We adopted a 
protocoP^ developed to address this issue, where the cell is subdivided into subvo- 
lumes (or compartments) and each reactant is treated as a different molecule in each 
subvolume, e.g., X^'^ 7^ X^^^ for the same compound X, where / and j refer to different 
subvolumes and the reaction propensities are altered accordingly. Each subvolume is 
assumed to be well-mixed - thus allowing for the use of the Gillespie algorithm. This 
protocol results in an increase in the number of reactions because (i) each reaction is 
treated independently in different subvolumes and (ii) additional "reactions" relating 
to the conversion of X^'^ to X^^^ (due to diffusion or translocation) are included, as 
described in detail by Bernstein^\ The nucleus (N), cytosol (C) and the mitochondrial 
membrane (M) are taken as three distinctive subvolumes here. No superscripts are 
appended for designating compounds in the cytosol, except for clarifying ambiguous 
cases. 

Statistical model checking. In order to test whether our stochastic model satisfies a 
given dynamical property with guaranteed confidence levels, we adopted a statistical 
model checking (SMC) based framework developed recently^^ SMC is a highly 
scalable simulation-based verification approach that can check system properties 
encoded as logic formulas. We use a strengthened bounded linear time temporal logic 
(BLTL) to specify qualitative or quantitative properties of interest. Specifically, let S be 
a finite set of real-valued variables and Tbe a positive integer. An atomic proposition 
(AP) in our logic is of the form x#y, where x and)/ are arithmetic expressions over real- 
valued variables in S, and #£{>,<,=,>,<}. The logic operators in our BLTL 
consist of A (and), V (or), (negation), and time bounded U (until), G (global), and F 
(future). The formulas of BLTL are defined as: (i) every AP as well as the constants 
true and false are BLTL formulas; (ii) if \J/, \J/' are BLTL formulas then -n/^ and \J/ V \J/' 
are BLTL formulas; (iii) if \J/, \J/' are BLTL formulas and t < Tis a positive integer then 
\J/U~^iJ/' and iJ/U*\J/' are BLTL formulas. The derived operators such as A, G~^ and 
are defined in the usual way. 

A trajectory in our model is a series of time- dependent states of the form a = (sq, 
h)y (5i> • •> meaning that the system jumps to state S/+i after staying in state for tf. 
We interpret the formulas of our logic at the finite set of time points T = {0,1,. . .,7j. A 
trajectory that ^satisfies a BLTL specified property (p at time f e T is written as a, t\ = 
(p. The semantics of the logic is defined as follows: (i) a, t\ = AP iff AP holds true in 
state St, (ii) and V are interpreted in the usually way; (iii) a, t\ = ^\J-^^' iff there 
exists /c' such that k" ^k,t= k' ^ Tand cr, f + = ijj'. (T,t + k"\ = \J/ for every 0 < k" 
< K\ (iv) (7, t\ = (//UV iii t + k ^ T and t + k\ = xjj'. a,t + k'\ = ij/ for every 
0 < ^' < A;. 

The statements we make are in the form of M| = Pr>o('A)> meaning that the 
probability that the system M satisfies a property ij/ is at least 9. For example, we 
express the property "caspase-3 level sustains once it reaches certain threshold" in 
Figure 3D as follows: 

Pr>o.95(C3<0.01nMU^^'^°HF^^^^(C3>0.3nM A G^4^^(C3>0.3nM)))) 



The BLTL expressions of other properties can be found in the Supplementary Text. 
Using SMC, the verification of such properties can be carried out approximately but 
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with guaranteed confidence levels and error bounds. Specifically, we verify each 
property using a sequential hypothesis test between the null hypothesis HO: p>6 + 5 
andthealternativehypothesisHl:p< 6 — (5, wherep is the probability of M satisfying 
\J/ and S specifies the indifference region. The strength of the test is determined by 
parameter a and which bound the type-I and type II errors, respectively. The test 
proceeds by generating a sequence of sample trajectories, Gi, (72,... A corresponding 
sequence of Bernoulli random variables Xi, %2> - • • are assumed, where = 1 if a, 0| = 
il/, otherwise Xk = 0. For each generated sample, we update the score a»„ by the 
following function: 

(6 - ^)ELi -'^'-^ (1 - (0 - (5))^"" ^"-1 "'^ 

««= 

(9 + ^)ELi -'•^'■-^ (1 - (0 + ^))('"- ELi 

where n is the number of generated samples. We accept hypothesis HO if co„ ^ 
(1 — jS)/a, and hypothesis HI if a»„ < jS/(l — a); otherwise, we draw another sample. 

Sensitivity analysis. Global sensitivity analysis is performed using a SMC-based 
multi-parametric sensitivity analysis (MPSA) method^^. We encode the steady state of 
caspase-3 (the model output) as a BLTL formula. The MPSA procedure involves 
drawing a representative set of samples from the parameter space. For each sampled 
combination of parameter values, we compute the objective value with respect to the 
BLTL property. The sampled parameter sets are classified into two classes using a 
threshold objective value. The sensitivities are then computed as the Kolmogorov- 
Smirnov statistics of cumulative frequency curves associated with the two classes. 
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